Carregando shapefile de são paulo:
library(rgdal)
## Carregando pacotes exigidos: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.0, released 2018/12/14
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ shared files: (autodetected)
## Linking to sp version:1.4-5
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
shp_sp <- readOGR(dsn = "data/shapefiles/SP_Municipios", layer = "SP", verbose = F)
plot(shp_sp)
Critério queen de vizinhança:
library(spdep)
## Carregando pacotes exigidos: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Carregando pacotes exigidos: sf
## Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
vizinhos_queen = poly2nb(pl = shp_sp, queen = TRUE,row.names = shp_sp@data$NM_MUN)
plot(shp_sp, border = "lightgray")
plot(vizinhos_queen,
coordinates(shp_sp),
add = TRUE,
col = "#33638DFF")
Matriz W Binária:
matrizW_queen <- nb2mat(neighbours = vizinhos_queen, style = "B", zero.policy = T)
colnames(matrizW_queen) <- shp_sp@data$NM_MUN
Gerando uma lista com várias matrizes W, cada uma com contiguidades de ordem 1 até 5:
vizinhos_queen_ordens <- nblag(neighbours = vizinhos_queen, maxlag = 5)
plot das vizinhanças com ordem 3:
plot(shp_sp, border = "lightgray")
plot(vizinhos_queen_ordens[[3]],
coordinates(shp_sp),
add = TRUE,
col = "#33638DFF")
Estabelecendo vizinhanças por contiguidade, critério rook (torre, não pega diagonais):
vizinhos_rook <- poly2nb(pl = shp_sp, queen = FALSE,row.names = shp_sp@data$NM_MUN)
plot(shp_sp, border = "lightgray")
plot(vizinhos_rook,
coordinates(shp_sp),
add = TRUE,
col = "#95D840FF")
Matriz W do critério rook:
matrizW_rook <- nb2mat(neighbours = vizinhos_rook, style = "B", zero.policy = TRUE)
colnames(matrizW_rook) <- shp_sp@data$NM_MUN
Shapefile do estado da Bahia:
shp_ba <- readOGR(dsn = "data/shapefiles/shapefile_ba/", layer = "ba_state", encoding = "UTF-8", use_iconv = TRUE, verbose=F)
vizinhanças por distâncias geográficas (as distâncias serão calculadas em km):
library(spdep)
vizinhos_distancias <- dnearneigh(coordinates(shp_ba), d1 = 0, d2 = 90, longlat = TRUE)
summary(vizinhos_distancias)
## Neighbour list object:
## Number of regions: 417
## Number of nonzero links: 11662
## Percentage nonzero weights: 6.706577
## Average number of links: 27.96643
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 5 6 6 2 4 4 7 10 6 16 6 10 7 10 4 12 8 17 18 20 15 11 12 11 7 5
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## 2 3 4 9 5 3 8 7 2 8 5 5 5 5 2 7 4 8 6 3 8 5 3 1 4 8
## 53 54 55 56 57 58 59 60 61 62 63 64 65 67
## 3 4 3 3 3 5 8 2 3 4 5 3 1 1
## 5 least connected regions:
## 113 393 408 409 415 with 1 link
## 1 most connected region:
## 240 with 67 links
Vizualizando as distâncias:
plot(shp_ba, border = "lightgray")
plot(vizinhos_distancias, coordinates(shp_ba), col = "#CC6A70FF", add =T)
Matriz W:
matrizW_distancias <- nb2mat(neighbours = vizinhos_distancias, style = "B")
colnames(matrizW_distancias) <- shp_ba@data$MUNICIPIO
rownames(matrizW_distancias) <- shp_ba@data$MUNICIPIO
library(rgdal)
shp_sc = readOGR("data/shapefiles/shapefile_sc", layer="sc_state", verbose=F)
Obrigando cada municipio ter exatamente 3 vizinhos:
lista_knear <- knearneigh(coordinates(shp_sc), longlat = TRUE, k = 3)
vizinhos_knear <- knn2nb(knn = lista_knear)
plot(shp_sc, border = "lightgray")
plot(vizinhos_knear, coordinates(shp_sc), add = T, col = "#13306DFF")
Matriz Binária W:
matrizW_knear <- nb2mat(neighbours = vizinhos_knear, style = "B")
colnames(matrizW_knear) <- shp_sc@data$NM_MUNICIP
rownames(matrizW_knear) <- shp_sc@data$NM_MUNICIP
matrizW_queen_linha <- nb2mat(vizinhos_queen, style = "W", zero.policy = T)
colnames(matrizW_queen_linha) <- rownames(matrizW_queen_linha)
São paulo tem 23 vizinhos no método queen:
sum(matrizW_queen["São Paulo",])
## [1] 23
Quando padronizamos por linha obrigamos a soma da linha dar 1:
sum(matrizW_queen_linha["São Paulo",])
## [1] 1
matrizW_queen_dupla_padr <- nb2mat(vizinhos_queen, style = "U", zero.policy = T)
colnames(matrizW_queen_dupla_padr) <- rownames(matrizW_queen_dupla_padr)
Quando padronizamos duplamente a soma por linha não da mais 1:
sum(matrizW_queen_dupla_padr["São Paulo",])
## [1] 0.006277293
Mas a soma da matriz inteira da um:
sum(matrizW_queen_dupla_padr)
## [1] 1
matrizW_queen_est_var <- nb2mat(vizinhos_queen, style = "S", zero.policy = T)
colnames(matrizW_queen_est_var) <- rownames(matrizW_queen_est_var)
Nesse caso a soma de matriz nos da a quantidade de municipios com vizinhos (ilha bela não tem vizinhos no método queen):
sum(matrizW_queen_est_var)
## [1] 644
Tipos:
Recuperando IDH de são Paulo e colocando no shapefile:
sp = read.csv("data/sp.csv")
shp_sp <- merge(x = shp_sp, y = sp, by.x = "CD_MUN", by.y = "codigo")
Para o cálculo da Estatística I de Moran, nosso algoritmo esperará como declaração um objeto de classe listw:
listw_queen <- mat2listw(matrizW_queen)
Estamos verificando a correlação de uma variável contra ela mesma defasada espacialmente (ou seja, comparando com os vizinhos) em relação a uma determinada observação:
Teste da Hipótese: o valor cálculado de I (0.2328220224) é estatisticamente ao valor esperado de I (0.00155521)? Sim, eles são estatisticamente diferente e rejeitamos a hipótese nula.
Valor esperado de I :
-1/(644-1)
## [1] -0.00155521
Estatística I de moran:
idhs = shp_sp@data$idh
idhs_scale = scale(idhs)
I = (t(idhs_scale) %*% matrizW_queen %*% idhs_scale /
t(idhs_scale) %*% idhs_scale) *
( (length(idhs)-1)/sum(matrizW_queen) )
I
## [,1]
## [1,] 0.232822
Como p-value = 6.695e-10 rejeitamos a hipótese nula e a estatsitica I de autocorrelação global é: 0.2328220224 - ou seja, existe autocorrelação global e ela é positiva.
moran.test(x = shp_sp@data$idh, listw = listw_queen, zero.policy = T)
##
## Moran I test under randomisation
##
## data: shp_sp@data$idh
## weights: listw_queen n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 10.092, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2328220224 -0.0015552100 0.0005393647
O Diagrama da Estatística I de Moran:
moran.plot(x = shp_sp@data$idh, listw = listw_queen, zero.policy = TRUE,
xlab = "IDH", ylab = "IDH Espacialmente Defasado")
Interpretação:
Diagrama de Moran com ggplot:
idh_defasada = lag.listw(x=listw_queen, var=shp_sp@data$idh,zero.policy = T)
df = data.frame(cidade = shp_sp@data$NM_MUN, idh=shp_sp@data$idh,idh_defasada)
df["idh_defasada_na_unha"] = matrizW_queen %*% shp_sp@data$idh
plotly::ggplotly(
df |>
ggplot(aes(label=cidade)) +
geom_point(aes(x=idh,y=idh_defasada), alpha=0.5,size=3) +
geom_smooth(aes(x=idh,y=idh_defasada), method="lm", se=F) +
geom_hline(yintercept = mean(df$idh_defasada),lty=2) +
geom_vline(xintercept = mean(df$idh),lty=2) +
annotate("text", x=0.66, y=18.5, label="Low-High") +
annotate("text", x=0.84, y=18.5, label="High-High") +
annotate("text", x=0.66, y=-1, label="Low-Low") +
annotate("text", x=0.84, y=-1, label="High-Low") +
theme_bw()
)
## `geom_smooth()` using formula 'y ~ x'
matrizW_queen_linha <- nb2mat(vizinhos_queen, style = "W",zero.policy = T)
colnames(matrizW_queen_linha) = rownames(matrizW_queen_linha)
listw_queen <- mat2listw(matrizW_queen_linha)
Estatística Moran Local com o uso da função localmoran:
moran_local <- localmoran(x = shp_sp@data$idh, listw = listw_queen, zero.policy = T)
idh_zscore = scale(shp_sp@data$idh)
moran_local_na_unha = rowSums(sweep(x=matrizW_queen_linha,MARGIN=2,STATS = idh_zscore, FUN = "*"))*idh_zscore
Juntando os resultados da Estatística Moran Local no dataset do objeto shp_sp:
moran_local_mapa <- cbind(shp_sp, moran_local)
library(gtools)
library(tmap)
quantile(moran_local_mapa@data$Ii, type = 5, probs = c(0,0.2,0.4,0.6,0.8))
## 0% 20% 40% 60% 80%
## -1.424372672 -0.123143140 0.002717804 0.121068801 0.458396014
moran_local_mapa@data <- moran_local_mapa@data %>%
mutate(faixa_quantis = factor(quantcut(x = Ii, q = 5)))
tm_shape(shp = moran_local_mapa) +
tm_fill(col = "faixa_quantis", palette = "-magma") +
tm_borders()
Considerando os quadrantes:
moran_local_mapa <- cbind(moran_local_mapa,
attr(x = moran_local, which = "quadr")[1])
tm_shape(shp = moran_local_mapa) +
tm_fill(col = "mean", palette = "-viridis") +
tm_borders(col = "gray")
Criando um vetor que contenha o centro das observações da variável idh ao redor de sua média:
quadrantes <- vector(mode = "numeric", length = nrow(moran_local))
centro_idh <- shp_sp@data$idh - mean(shp_sp@data$idh)
centro_moran_local <- moran_local[,1] - mean(moran_local[,1])
Vamos manter na base somente os I locais significantes. Enquadrando nossas observações em seus respectivos quadrantes:
sig <- 0.05
quadrantes[centro_idh > 0 & centro_moran_local > 0] <- "HH"
quadrantes[centro_idh > 0 & centro_moran_local < 0] <- "HL"
quadrantes[centro_idh < 0 & centro_moran_local > 0] <- "LH"
quadrantes[centro_idh < 0 & centro_moran_local < 0] <- "LL"
Ajustando a presença da observação em razão de sua significância estatística:
quadrantes[moran_local[,5] > sig] <- "Estatisticamente_não_significante"
moran_local_mapa@data["quadrantes"] <- factor(quadrantes)
Plotando os quadrantes de forma espacial (versão ‘default’):
tm_shape(shp = moran_local_mapa) +
tm_polygons(col = "quadrantes",
pal = c(HH = "darkred",
HL = "red",
LH = "lightblue",
LL = "darkblue",
Estatisticamente_não_significante = "white")) +
tm_borders()
## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
Recarregando dados da Bahia:
library(rgdal)
library(tidyverse)
library(spdep)
shp_ba <- readOGR(dsn = "data/shapefiles/shapefile_ba/", layer = "ba_state",
encoding = "UTF-8", use_iconv = TRUE, verbose=F)
vizinhos_distancias <- dnearneigh(coordinates(shp_ba), d1 = 0, d2 = 90, longlat = TRUE)
matrizW_distancias <- nb2mat(neighbours = vizinhos_distancias, style = "B")
colnames(matrizW_distancias) <- shp_ba@data$MUNICIPIO
rownames(matrizW_distancias) <- shp_ba@data$MUNICIPIO
idh_bahia = read.csv("data/idh_bahia.csv")
Adicionando dados de IDH da bahia ao shapefile:
idh_bahia$Codigo = as.character(idh_bahia$Codigo)
shp_ba@data = shp_ba@data |> left_join(idh_bahia, by = "Codigo")
Transformação do objeto matrizW_distancias em um objeto de classe listw:
listw_dist <- mat2listw(x = matrizW_distancias)
Calculando a Estatística G de Getis e Ord:
g_local <- localG(x = shp_ba@data$idh, listw = listw_dist)
Juntando as informações do objeto g_local ao nosso shapefile:
mapa_G <- cbind(shp_ba, as.matrix(g_local))
mapa_G@data = mapa_G@data |> rename(estistica_g = 4)
Plotando a Estatística G de forma espacial:
library(tmap)
tm_shape(mapa_G) +
tm_fill("estistica_g", palette = "-RdBu") +
tm_borders()
Estradificação por 8 e melhorado para daltônicos:
library(gtools)
mapa_G@data = mapa_G@data |> mutate(faixa_quantis = factor(quantcut(x = estistica_g, q = 8)))
tm_shape(mapa_G) +
tm_fill("faixa_quantis", palette = "plasma") +
tm_borders()